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ABSTRACT 


Supersonic nozzles which operate at low Reynolds 
numbers and have large expansion ratios have very thick 
boundary layers at their exit. This leads to a very strong 
viscous/inviscid interaction upon the flow within the noz- 
zle and the traditional nozzle design techniques which 
correct the inviscid core with a boundary layer displace- 
ment do not accurately predict the nozzle exit conditions. 
In addition, if the nozzle exit density becomes low enough 
rarefaction effects in the form of velocity slip and tem- 
perature jump at the wall must be accounted for. 

The present work vised a full Navier-Stokes code 
(PARC 2D) to compute the nozzle flow field. Grids were gen- 
erated using the interactive grid generator code TBGG . All 
computations were made on the NASA MSFC CRAY X-MP computer. 
Comparison was made between the computatiohs and in— house 
wall pressure measurements for CO 2 flow through a conical 
nozzle having an area ratio of 40. Satisfactory agreement 
existed between the computations- -and measurements for a 
stagnation pressure of 29.4 psia and stagnation temperature 
of 1060 °R. However, agreement did not exist at a stagna- 
tion pressure of 7.4 psia. Several reasons for the lack of 
agreement are possible. The computational code assumed a 
constant gas gamma whereas gamma for CO 2 varied from 1.22 
in the plenum chamber to 1.38 at the nozzle exit. The com- 
putations were performed assuming adiabatic, no-slip walls, 
both of which may not be correct. Finally, it is possible 
that condensation occurred during the expansion at the 
lower stagnation pressure. The next phase of the work 
will incorporate variable gamma and slip wall boundary 
conditions in the computational code. 
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INTRODUCTION 


Low Reynolds number, large expansion ratio, and super- 
sonic nozzles are frequently used for control of space- 
craft. An accurate knowledge of their thrust is required 
for designing the spacecraft control system. Recently 
the contamination to the spacecraft from the nozzle 
exhaust has become a concern and accurate computations of 
the exhaust plume have been required. The flow in such 
nozzles possesses strong viscous/inviscid interactions at 
their exit due to the thick boundary layers. Traditional 
nozzle design techniques, such as the use of the method of 
characteristics to calculate the inviscid core and bound- 
ary layer theory to compute the displacement thickness, 
fail to predict the strength of the viscous/inviscid 
interaction. Therefore, it is necessary to use a full 
Navier-Stokes code for this purpose. 

The present work used an existing code that incorpo- 
rates the full Navier-Stokes equations and viscous stress 
terms to calculate the flow fields within a conical nozzle 
having an area ratio of 40. Comparison was made with in- 
house measurements that had been performed on this nozzle 
using CO 2 as the gas. Tests were performed at several 
stagnation pressure levels, the lowest resulting in a exit 
wall Knudsen number of 0.06. Under these conditions 
slip flow conditions could be expected to exist at the 
nozzle exit. This presents an additional complication to 
the computation of the flow field. The PARC2D code, 
developed by Sverdrup Technology, Inc., AEDC Group (ref. 6) 
has demonstrated the ability to calculate such nozzle flow 
fields and their plumes. The slip flow could be included 
as a new boundary condition subroutine within the program. 
Therefore, it was chosen to perform this task. 
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OBJECTIVES 


The objectives of this work were to: 

1) generate a computational grid using the interactive 
program TBGG on the MSFC VAX 785 computer, 

2) modify the output file from TBGG and generate 
additional input files for PARC2D, copy them to a tape 
and transfer the files to the EADS system, 

3) install and make operational on the MSFC CRAY machine 
the Navier-Stokes code PARC2D, 

4) . couple the grid generator output file to the PARC 2D 
code, 

5) run PARC2D for conditions that matched those of in- 
house measurements of CO 2 flow through a 40:1 area ratio 
nozzle 

a) determine h<aw to converge the program solution, 

b) compare' the computed results with the measure- 
ments as the value of gamma was varied, 

c) compare the computed results and measurements as 
the grid and downstream boundary condition were 
varied , 

6) learn how to modify PARC2D to include slip wall bound- 
ary conditions, 

7) learn how to run PARC 2D to obtain the nozzle plume 
flow field and to calculate the force on a nozzle end 
plate. 
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GRID GENERATION 


The first step to calculating the flow in a nozzle is 
to generate an appropriate computational grid. Two grid 
generator programs acquired from Arnold Engineering Develop- 
ment Center were placed on the MSFC CRAY machine but some 
of the routines were found to be incompatible with the IBM 
front end and, rather than rewrite these portions of the 
programs, it was decided to use a grid generator that had 
previously been installed on the VAX system. 

The grid generator that was used is called TBGG and 
was developed by Smith and Wiese at Langley Research Center 
(ref. 1) . It is interactive in the 4010 mode and can be 
used to generate body-fitted grids algebraically. The 
number of grid points can be changed, with the maximum 
possible number controlled by a parameter IPX. This 
number was initially 100 but was later changed to 200. 

The program starts from a RESTART file that had been pre- 
viously created (for given geometry) or from a file DATANEW 
that is created and contains a listing of the boundary 
coordinates. The grid lines can be concentrated at differ- 
ent regions of the flow field using several functions. 

For the nozzle flow fields the bi-exponential function was 
used for the top and bottom surfaces and the connecting 
curve. TBGG generates two output files, RESTART, which is 
an unformatted file that can be used to start the program 
next time, and GRIDOUT, which is a formatted output file. 

A modification of this latter file was used as part of the 
necessary input to the Navier-Stokes code. 

Computations were performed on a nozzle configuration 
that was used for MSFC in-house measurements, using C0 2 as 
the gas. The major geometric parameters for the nozzle and 
upstream plenum chamber are shown in Figure 1. This geome- 
try was introduced into TBGG as a discrete number of 
dimensional points. The bi-exponential function was used 
to cluster the grids perpendicular to the nozzle axis 
(x-axis in Figure 2) in the region of the nozzle throat by 
choosing Ki = 0.57 and K 2 = -2.0 for both top and bottom 
curves (see ref. 1). For the connecting curve Ki = 0.50 
and K 2 = 2.0 were chosen to place the grid points near the 
nozzle walls and symmetrically about the centerline 
(Figure 3) . Because the Navier-Stokes code only required 
grid points from the centerline to the wall, 99 grid points 
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Figure 1. 


Geometry of conical nozzle and plenum 
chamber. All dimensions are in inches. 


A e /A* =40. 


were chosen in the y-direction to place one on the center- 
line and 100 were chosen in the x-direction. Enlargements 
of the grid next to the wall in the throat region and at 
the nozzle exit are shown in Figures 4 and 5, respectively. 

The output file GRIDOUT gave the arrays 
[ (X (J, K) , K = 1, KMAX) , J = 1, JMAX] and the same for Y(J, K) , 
where J is the grid index in the x-direction (JMAX = 100) 
and K is the index in the y-direction (KMAX = 99) . As 
stated above, the Navier-Stokes code required grid points 
only from the nozzle centerline to one wall and it required 
that they be ordered in the opposite fashion. Thus, the 
code required [ (X (J, K) , J = 1, JMAX) , K = 1, KMAX] where now 
KMAX was 50. A program was written to accomplish this plus 
generate the other necessary input arrays (see next section) . 
This file was copied to magnetic tape and read onto the 
EADS system. 

The grid just described was discovered to be inadequate 
(see next section) . Therefore, a second grid was gener- 
ated using the geometry shown in Figure 6. This geometry 
contained one inch less plenum chamber, which was added as 
an exhaust chamber downstream of the nozzle exit. It was 
also decided to increase the number of grid points to 200 x 
199. This required changing the parameter statement, as 
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Figure 2. 


Schematic of geometry used for 100 x 99 grid 
All (x,y) coordinates are defined with 
respect to this geometry. 



Figure 3. Portion of 20 x 20 grid applied to geometry 
above, illustrating the concentration of 
grid lines at outer walls and nozzle throat. 



Figure 5. 


Blow-up of 100 x 99 grid next to nozzle wall 
at exit plane. 
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Figure 6. 


Schematic of geometry used for 200 x 199 grid, 
with nozzle exiting into a tube. 



Figure 7. 20 x 20 grid applied to geometry above, 

illustrating concentration of grid lines at 
the outer walls and nozzle throat. 
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stated earlier. Again the grids lines were concentrated 
near the nozzle throat with Ki = 0.33, K 2 = -2.0 and near 
the walls with = 0.50, = 2.0 (Figure 7). Problems 

arose when trying to get a useful grid that would flow 
from the nozzle into the exhaust chamber. A perpendicular 
end plate failed to give continuity to the grid lines 
(Figure 8) and since it was desirable to have one set of 
grid lines perpendicular to the nozzle axis (to be able to 
plot properties across a nozzle section) a 1/16 inch 
radius connecting the nozzle exit with a perpendicular end 
plate was also unsatisfactory (Figure 9). However, adding 
a slope to the end plate did yield a useful grid (Figure 10) 
This grid, modified to 200 x 100, was the one that was used 
for the latter calculations. 


NAVIER- STOKES COMPUTATIONAL CODE 


The Navier-Stokes code, PARC 2D, is a modification of 
the ARC2D code that was developed by Pulliam and Steger at 
NASA Ames. Space does not permit more than a cursory 
description of that code but more details are given in 
references 2 to 4. It uses a thin layer approximation to 
the Navier-Stokes equations with parabolized viscous stress 
terms. It is written in strong conservative form in cur- 
vilinear coordinates, utilizing the delta form (ref. 7). 

The noniterative implicit approximate factorization scheme 
of Beam and Warming (ref. 5) is used with fourth order 
artificial dissipation. 

The PARC2D code is a modification of the ARC2D code by 
Sverdrup Technology, Inc., AEDC Group which removed the thin 
layer approximation and the approximate stress terms (ref. 6) 
It is fully elliptic, requiring closed boundary conditions 
and an initial condition everywhere in the flow field. The 
code is modular and fully vectorized. It assumes that the 
gas is a perfect gas with constant gamma and Prandtl number 
but uses the Sutherland viscosity law for the temperature 
variation of viscosity. An initial test over 400 iterations 
indicated that the vectorization decreased the CPU time on 
the MSFC EADS system by a factor of 5.22. 

Some of the important required input quantities to the 
code are listed in Table I. The Reynolds number is based 
upon a reference sound speed and length. For the present 
calculations the stagnation conditions were taken as refer- 
ence and the reference length was one inch, since the grid 
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TABLE I 


OPERATIONAL ASPECTS OF PARC 2D 


Nondimensionalization: q = q/a r ; p = p/y p r ; p = p/p r 



Re 


p r a r L 


Hr 


q = speed, a = sound speed, 
( ) = reference 


Namelist inputs: y, Re, Prandtl number, P r , T r , 

J max , K max selection from adia- 
batic/constant wall temperature, 
viscous/inviscid, axisymmetric/ 
two-dimensional, laminar/ 
turbulent 


Initial condition input: Iteration number, gamma 

arrays X(J, K) , Y(J, K), p(J,K), 
p u(J, K) , p V ( J , K) , E ( J, K) 
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coordinates were calculated in inches. Although several 
diff eren t conditions were used as initial conditions it was 

found that th© fa^stost conv©3rg©nc© occurred by taking u = 

v = o and p and E equal to their stagnation values, 
f*®*f P = ^ • E = 1/y(Y — 1). These values were assigned 
initially throughout the flow field. Other options used are 
as follows: adiabatic wall, axisymmetric, viscous, laminar. 
Free boundaries (allowing inflow or outflow) were assumed at 
the upstream and downstream boundaries in Figures 2 and 6 , 
and the other boundaries were axis of symmetry and no slip/ 
adiabatic wall. it would have been possible to have taken 
the upper and lower portions of the chamber downstream of 
the nozzle in Figure 6 as outflow boundaries also but that 
was not done. A parameter statement for NX, NY and NM had 
to be changed to allow JMAX = 200. Note that NM must be at 
least as large as the largest of the other two parameters. 

Boundary conditions must be imposed upon all of the 
boundaries. For the free boundaries this consists of a 
specification of the pressure. The stagnation pressure was 
imposed upstream but the downstream boundary condition posed 
a problem. For flows with strong viscous/inviscid inter- 
action the pressure cannot be expected to be constant across 
even a contoured nozzle and even greater pressure variations 
would exit across the exit of a conical nozzle. However, a 
constant downstream pressure must be specified. It is 
usually prudent to specify a pressure somewhat lower than 
the minimum expected downstream pressure. No boundary con- 
ditions are required on the axis nor on the wall. If constant 
temperature wall conditions are assumed then the temperature 

be specified. Note that the entire wall does not have 
to be assumed to be at the same temperature. 

Placing the downstream boundary conditions at the nozzle 
sxit plane produced unacceptable stagnation pressure fluctua- 
tions on the centerline upstream at the exit (Figure 11) . 

It was for this reason that the second geometry (Figure 6) 
was used for the latter calculations. The boundary condition, 
however, had only a small effect upon the other results, such 
as the exit Mach number, as can be seen from the two calcula- 
tions given in Table II for y = 1.33 and Po = 29.4 psia. 

These two Mach numbers differ by only 0.3 percent. 

An efficient convergence procedure for these nozzle 
problems was discovered. The steps are as follows: 

1) Initially set q = 0, p = 1 and E = 1/y(y~1) every- 
where in flow field. 
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Figure 11. Centerline stagnation pressure approaching the 
exit plane. p Q = 29.4 psia, T q = 1060 °R, 

Y = 1.3, Re = 68,115, 100 x 50 grid. 


2) Set initial parameters as follows: DIS2 = 0.2, 

DIS4 = 0.35, PCQMAX = 10.0, DTCAP =5.0. Run until 
the axial velocity is positive everywhere in the 
plenum chamber. 

3) Slowly reduce DIS4 to 0.25 (all other parameters 
constant) . 

4) Slowly reduce DIS2 to 0.00. 

5) Check DT and set DTCAP to about one-half of the 
minimum DT for the last series of interations. 

Then run until L2 reaches an acceptable value 

(10 -8 to 10" 9 ) . 

In this discussion DIS4 and DIS2 are parameters related to 
the fourth order and second order dissipation, respectively, 
PCQMAX sets the maximum change in any variable during an 
iteration, DT is the time step and DTCAP is the maximum 
allowable time step. L2 is a convergence measure. 
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Convergence is also improved by minimizing the amount 
of plenum chamber that is included in the computational 
domain. Small disturbances in this low velocity region 
damp very slowly. 


RESULTS 


Computations were made using the test conditions that 
corresponded to the MSFC in-house measurements on the noz- 
zle described in Figure 1. The test gas was C0 2 , stagna- 
tion temperature was 1060 °R and the stagnation pressure 
had two values , 29.4 and 7.4 psia. Since only wall pressure 
measurements were made, that was the only parameter that 
could be used for comparison with the computations. The 
computational code assumed a constant gamma whereas the 
gamma for C0 2 varied from about 1.22 in the plenum chamber 
to 1.38 at the exit. Therefore, exact comparison could 
never be expected to occur between the measurements and 
computations . 

A summary of some of the results is given in Table II. 
Many of the initial computations were performed at an 
erroneous Reynolds number but their results are included 
for completeness. Comparison of computed and measured wall 
pressures is given in Figure 12. The agreement at the 
higher stagnation pressure is very good (also see Figure 13 
which is drawn to magnify the differences) and could be 
improved by increasing the assumed value of y. 

The effect of y on the results is shown by the two 
computed results at the lower pressure. The y = 1.33 
computations compare more favorably with the measurements 
than do those at y = 1.30 but the agreement is still not 
good. Several reasons are possible for the lack of agree- 
ment, in addition to the need for a variable y computational 
code. These include the possibility that the adiabatic wall 
boundary condition is not applicable and computations at a 
constant wall temperature should be performed to examine 
that possibility. In addition, there is a strong possibility 
that condensation of the C0 2 was occurring during the expan- 
sion and agreement would not be expected until the stagna- 
tion temperature was raised to eliminate all possibility for 
condensation. Finally, at the exit the Knudsen number is 
about 0.06, based on a mean free path at wall conditions 
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and a crude estimate of the displacement thickness. Slip 
is expected to occur under such circumstances and the slip 
wall boundary condition formulated by the author (ref. 8) 
should be implemented rather than the no-slip condition. 

Typical profiles of Mach number and static temperature 
across the nozzle at a location near the nozzle exit are 
shown in Figure 14. The large temperature gradients near 
the wall are caused by the assumed adiabatic wall condition. 

Very large wall pressure gradients were observed at 
the nozzle exit. For the original grid. Figure 2, this 
influence extended upstream several grid points in the sub- 
sonic boundary layer. This was caused by the requirement 
of the code to exactly match the prescribed downstream 
boundary pressure at the wall and was one of the reasons 
for going to the configuration shown in Figure 6. However, 
a rapid expansion still existed right at the exit because of 
the details of the grid that was generated (Figure 10) . 

In this case the expansion only occurred over the last grid 
point. This problem can be eliminated by slightly extend- 
ing the nozzle and using the computed properties only up 
to the actual nozzle exit plane. 

The usefulness of the computational results can be 
greatly enhanced by the implementation of plotting 
facilities that are part of the EADS system and then Mach 
number contours, for example, can be drawn for the entire 
nozzle. 
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TABLE II 


SUMMARY OF RESULTS 
(Adiabatic wall) 

Convergence information: p Q = 29.4 psia 

Number 

Y Re Grid Iterations L^ CPU Time 

1.4 70/686 100 x 50 23,000 4.5 x 10~ 9 40 min 

1.3 664,341 200 x 100 9,500 3.3 x 10~ 9 64 min 


Exit properties at centerline: 


P Q (psia) 

Y 

Re 

Grid 

Centerline 
Mach Number 

29.4 

1.30 

68,115 

100 x 50 

4.295 


1.33 

68,896 


4.409 


1.40 

70,686 


4.724 


1.33 

68,896 

200 x 100 

4.395 

• 

1.31 

666,891 


4.652 

7.4 

1.40 

17,792 

100 x 50 

4.188 


1.30 

168,496 

200 x 100 

4.454 


1.33 

170,429 


4.601 


1.40 

174,856 


4.967 


1.33 

170,429 


4.671 
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Figure 12. Comparison of measured and computed nozzle 
wall pressure distributions. Computations 
used 200 x 100 grid. T Q = 1060 °R. 
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Figure 13. Comparison of measured (CO2) and 
nozzle wall pressures. p Q = 29. 
T 0 = 1060 °R. Calculations used 
and 200 x 100 grid. 


computed 
psia , 

Y = 1.31 
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y (inches) 


Figure 14. Computed Mach number and static temperature 
profiles across the nozzle at x = 4.342 
inches. p Q = 29.4 psia, T 0 = 1060 °R, 

Y = 1.31, Re = 666,891, adiabatic wall. 
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CONCLUSIONS AND RECOMMENDATIONS 


The following conclusions can be made from results of 
this study. 

1) The grid generator code TBGG provided a good means 
generating useful computational grids for the 
Navier-Stokes computations of the nozzle flow. 
However, it had several drawbacks. First, it was 
inconvenient to have to generate the grid on the 
VAX, transfer it to a tape and place it on the 
CRAY. Second, it was difficult to concentrate the 
grid points in all of the regions of interest. 
Finally, it would be difficult to use this grid 
generator for more complicated geometries and grid 
generators that can patch various regions together 
are available for such circumstances. 

2) The PARC 2D code yields accurate solutions for the 
flow field in supersonic nozzles at low Reynolds 
numbers where large viscous/ inviscid interaction 
exists. Because the code is modular it can be 
easily modified. An efficient convergence pro- 
cedure was developed for these nozzle flows. The 
amount of plenum chamber included in the computa- 
tion should be minimized and an exhaust chamber 
should be provided downstream of the nozzle exit 
to increase the accuracy and minimize the number 
of iterations. 

A number of recommendations for additional work or 
improvements to the code can be made. 

1) The constant gamma restriction to the PARC 2D code 
should be removed. Sverdrup Technology, Inc., 

AEDC Group, is presently modifying the code to 
accommodate variable gamma. 

2) Additional computations should be made at other 
values of gamma, especially for the low stagnation 
pressure test, to see if improved agreement of the 
computations with the measurements could not be 
achieved. Computations should also be performed 
at constant wall temperature rather than an 
adiabatic wall to see the influence of that 
boundary condition upon the flow. 
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3) Subroutines should be written to calculate bound- 
ary layer thicknesses, etc., to allow comparison 
with other computational techniques. 

4) The plotting routines available on the CRAY 
system should be utilized to better understand 
the nature of the solutions that have been 
obtained. 

5) Objectives numbers 6) and 7) were not accomplished 
due to lack of time. These should be implemented, 
especially the slip boundary condition. That con- 
dition will be required, along with the variable 
gamma, before the computations will satisfactorily 
agree with the measurements at the lowest stagna- 
tion pressure. Also, flow should be allowed to 
leave the sides of the downstream region to more 
accurately model the nozzle plume. 
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